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I.  INTRODUCTION 


A . Background 

The  Northrop  Corporation,  under  Air  Force  funding,  has  deveJoped 
a finite  element  digital  computer  code,  called  BR-1,  for  predicting  the 
inelastic,  large  deflection,  transient  response  of  combat  aircraft  skin- 
rib-stringer  structures  when  subjected  to  internal  air  blast  loading. 

The  finite  elements  considered  are  flat  rectangular  plates  and  beam 
stiffeners.  The  theory,  user's  manual  and  code  listing  are  given  in 
References  1 and  2.  The  Air  Force  Flight  Dynamics  Laboratory  wanted 
the  BR-1  code  modified  so  that  it  could  be  used  to  predict  the  response 
of  aircraft  fuel  tank  walls  when  subjected  to  fluid  pressures  due  to 
projectiles  passing  through  the  fuel  in  the  tank.  The  intense  pressure 
and  momentum  in  the  fuel  due  to  the  penetrating  projectile  is  referred 
to  as  the  hydraulic  ram  loading.  This  report  describes  the  modifica- 
tions to  the  IBM  version  of  the  BR-1  code  to  account  for  the  fluid 
(fuel)  - structure  (tank  wall)  interaction  that  occurs  when  bullets  and 
metal  fragments  penetrate  into  aircraft  fuel  tanks.  The  modified  code 
is  called  BR-1HK , The  interaction  between  the  compressible  fluid  and 
the  structure  is  approximated  by  the  piston  theory.  The  code  can  also 
be  used  for  many  other  compressible  fluid-structure  interaction  problems. 

B.  Piston  Theory 

The  total  nonlinear  problem  of  the  response  of  a tank  containing  a 
fluid  and  subjected  to  a high  speed  penetrating  projectile  is  extremely 
complex  and  presently  defies  analytical  treatment.  In  general,  the 
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equations  for  the  fluid  stresses  and  motion  are  coupled  to  those  for  the 
wall  stresses  and  motion  due  to  the  continuity  at  the  fluid-structure 
interface  (3).  One  procedure  for  approximating  the  fluid-structure  in- 
teraction phenomenon  is  the  piston  theory  (4).  This  theory  has  been  in 
use  since  the  early  1940' s when  it  was  applied  to  the  study  of  the 
effect  of  underwater  explosions  on  ship  plates.  It  provides  the  correct 
solution  to  the  one -dimensional  propagation  of  stresses  in  an  acoustic 
medium  due  to  a moving  boundary.  Several  recent  studies  have  been  made 
to  determine  its  accuracy  when  applied  to  two  dimensional  fluid-structure 
interaction  problems  (4,5). 

Application  of  the  piston  theory  to  the  interaction  problem  allows 
the  structure  equations  and  fluid  equations  to  be  uncoupled.  The  response 
of  the  wall  is  computed  using  the  conventional  structural  response  equa- 
tions, with  the  normal  pressure  on  the  wall  p given  by 

p = Po  + DC  (v±  - w)  (1) 

where  pq  and  v^  are  the  incident  pressure  and  velocity  of  the  fluid  at 

the  wall  respectively,  o is  the  fluid  density,  c is  the  acoustic  velocity 

in  the  fluid,  and  w is  the  wall  velocity*.  The  pressure,  pq,  and 

velocity,  v^,  are  the  pressure  and  velocity  that  would  exist  in  the 

fluid  if  the  interface  was  not  there,  i.e.,  p and  v.  do  not  contain 

’ o i 

any  "local"  reflected  effects.  However,  effects  on  p and  v.  due  to 

o i 

earlier  reflections  from  other  walls  and  free  surfaces  should  be  con- 
sidered. In  other  words,  pQ  and  v^  are  the  loading  components  due  to 

*A  dot  above  a variable  denotes  a derivative  with  respect  to  time. 
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the  free  field  and  the  scattered  effects.  The  loading  component  due  to 

the  wall  velocity  w is  called  the  radiation  pressure. 

C.  The  NWC  Hydraulic  Ram  Computer  Code 

In  order  to  use  the  piston  theory  to  compute  the  tank  wall  response, 

it  is  necessary  to  know  the  incident  fluid  pressure  p and  velocity  v 

o i 

over  the  entire  fluid -wall  interface  as  a function  of  time.  In  conjunc- 
tion with  this  project  Lundstrom,  at  the  Naval  Weapons  Center,  has 
developed  a digital  co.nputer  code  that  predicts  the  fluid  pressures  and 
velocities  pQ  and  v throughout  a rectangular  body  of  fluid  due  to  a 
penetrating  ballistic  projectile.  The  model  is  based  upon  replacing 
the  projectile  by  a line  of  sources  whose  strength  is  determined  by  an 
energy  balance  between  the  kinetic  and  potential  energy  of  the  fluid  and 
the  energy  loss  due  to  drag  forces  on  the  projectile.  Reflections  from 
the  structure -fluid  interface  are  accounted  for  by  considering  the  fluid 
boundary  to  be  either  stress  free  or  rigid*.  An  extensive  series  of 
tests  were  performed  at  the  Naval  Weapons  Center  to  obtain  detailed 
pressure  measurements  for  a variety  of  projectiles  under  a wide  range  of 
impact  conditions.  This  data  allowed  the  selection  of  important  para- 
meters such  as  tumbling  distance,  Jacket  stripping,  etc.,  to  be  entered 
into  the  code.  A description  of  the  code  and  the  instructions  for 
operation  are  given  in  Reference  7*  This  code  provides  the  values  for 
pQ  and  v^  at  user  specified  locations  over  the  structure -fluid  interface 
for  the  time  span  of  interest. 

* A study  of  the  one-dimensional  reflection  of  step  pressure  waves  from 
typical  aircraft  fuel  tank  walls  indicates  that  the  stress  free  surface 
provides  the  more  accurate  approximation  (6) 
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II.  MODIFICATION  OF  THE  BR-1  CODE 


A.  Incorporation  of  the  Piston  Theory 

The  BR-1  code  has  an  option  for  the  user  to  input  a time  varying 
pressure  on  each  panel  element.  In  the  piston  theory  this  pressure  is 
the  pQ  + pcv^  of  Eq.  1.  The  other  contributor  to  the  wall  loading  given 
by  Eq.  1 is  w,  the  wall  velocity.  Since  the  BR-1  code  does  not  include 
damping  effects,  it  is  necessary  to  add  the  damping  term  pew  to  the 
equations  of  motion. 

The  BR-1  code  solves  the  set  of  equations  (Eqs.  1 and  2,  Ref.  1) 

[M]  {q*J  = [P)  - [P)  - [H]  (q.)  - {CJ  (2) 


for  the  vector  of  global  nodal  generalized  displacements  {q*}  as  a function 
of  time.  These  generalized  displacements  define  the  motion  of  the  walls. 
The  vector  {F)  consists  of  global  generalized  external  and  body  forces 
at  the  nodes  of  the  elements.  The  matrix  [M]  is  the  mass  matrix. 

The  wall  pressure  p given  by  Eq.  1 causes  external  forces  at  the 
nodes.  The  external  generalized  forces  at  the  nodes  of  each  element  in 
the  local  coordinate  system,  {f},  is  given  by  (Eq.  A-l+7,  Ref.  1) 


tf)  ■ 


clA 


(3) 


where  [nb]  is  the  transpose  of  the  matrix  of  shape  functions  [N], 
evaluated  at  the  surface  of  the  element,  Aout  is  the  surface  of  the  ele- 
ment, and  {Tq)  is  the  vector  of  applied  surface  tractions  and  moments. 
The  order  of  {T  } is  a 5x1  vector.  Due  to  the  fluid  pressure  loading 
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i*u 


0 
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where  the  subscripts  denote  the  coordinates,  (T^  is  normal  to  the  ele- 
ment), and  p is  given  by  Eq.  1.  The  fourth  and  fifth  elements  of  {Tq} 
correspond  to  applied  moments  per  unit  middle  surface  area  of  the 
element,  and  are  zero  here.  The  numerical  intergration  of  Eq.  3 can 
be  accomplished  by  Gaussian  quadrature.  However,  the  {f}  is  obtained 
in  the  BR-1  code  in  a more  approximate  way  by  using  a lumping  approach 
at  the  nodes  ci”  the  element,  as  is  done  for  the  mass  matrix  evaluation, 
in  order  to  save  computation  time.  Thus,  according  to  Eq.  B-92  of 
Ref.  1,  the  external  force  vector  at  the  rth  node  of  the  nth  element 
is  given  by 

®2 
-0i 

i 
o 
0 

where  pn  is  the  magnitude  of  the  pressure  on  the  element*,  and  Dnr  = 

V(1+8H>  nr  where  and  0g  are  the  fourth  and  fifth  local  general- 
ized displacements  at  the  rth  node.  They  appear  here  because  the 
pressure  is  defined  in  BR-1  as  the  pressure  normal  to  the  deformed 
surface.  The  quantities  (xg-x^)  and  (y^-y2)  are  the  dimensions  of  the 
rectangular  element. 

The  pressure  pR  in  the  piston  theory  is  given  by  Eq.  1,  i.e. 

p = p + (pc)  v.  - (pc)  w (6) 

rn  on  n in  y n nr  ' 7 


* The  assumption  is  made  in  the  programming  of  BR-1  that  the  pressure 
is  uniform  over  each  element.  This  is  contrary  to  the  theoretical 
presentation  where  the  pressure  is  defined  at  each  node  point. 


where 


(pc) 


pc  if  the  nth  element  is  in  contact  with  the  fluid 


n 


0 if  the  nth  element  is  not  in  contact  with  the  fluid. 

The  variables  p and  v.  can  be  determined  from  the  NWC  computer  code 
on  m 

for  each  element  for  the  time  span  of  interest  prior  to  the  computation 

of  the  wall  response.  This  data  is  then  input  as  the  known  external 

pressure.  The  variable  w is  an  unknown  dependent  variable  and  is  part 

nr 

of  {q*}.  hence,  it  must  be  incorporated  into  the  equations  of  motion, 
Eq.  2. 

The  ff  1 due  to  w is  given  by 
1 nrJ  nr  ° J 


{f  1 
1 nr- 


•M  Mn 


(pC)n  Wnr 


nr 


1 

0 

0 


(7a) 


nr 


The  global  force  vector  at  the  rth  node,  fF]r,  is  related  to  {fnr}  in 


the  form 


{FI  = E [J  ] ff  } 

1 Jr  ft  L nJ  1 n«J 


Of — *r 


(7b) 


where  [J  ] is  the  transformation  matrix  from  the  global  coordinate 
system  to  the  nth  element  local  coordinate  system,  g means  a summation 
over  all  elements  containing  the  node  r,  and  a — r means  node  <*  corre- 
sponds physically  to  node  r.  Since  the  .wall  velocity  in  Eq.  7a  is  given 
in  terms  of  the  local  coordinates,  it  must  be  converted  to  global 
coordinates.  Thus,  according  toEqs.  A-77  ana  B-3  of  Reference  1, 


w = I J j fq*} 
nr  L nJ 


(8) 


3 

where  |_J^J  1®  the  third  row  of  [J  ].  Thus,  the  global  external  force 
vector  at  the  rth  node  (F}r  becomes 


Note  that  {F}^  is  nonlinear  since  9^  and  0^  are  Part  of  If  the 

rotations  ©1  and  are  neglected  in  the  computation  of  {F}r,  i.e.,  if 
the  pressure  is  not  truly  normal  to  the  deformed  surface,  {F}  for  the 
total  R nodes  of  the  structure  can  be  given  in  the  form 


where 


{F}  * - [D]  {4*} 


CD]  = 


and  [D]r  is  a 6x6  matrix  given  by 


(10a) 


(10b) 


(10c) 
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B.  Method  of  Solution 


The  BR-1  code  solves  for  {q*}  at  discrete  points  in  time  using  the 
explicit  finite  difference  scheme  (Eq.  A-109,  Ref.  l) 


(Aq*}.  = {Aq*}  + At  {q*}  + (At)  (q*}. 

l+l  i 1 xi 


(11) 


where  At  is  the  time  interval  between  two  time  points,  i.e. 


and 


* = Vi  - 


= (q*J  . (q.) 

i+1  ^i+l  1 


The  acceleration  fq*},  is  obtained  from  Eq.  (2)  in  the  form 

ti 


(12) 


fq*}t  = [M]-1  CC]  (13) 

Li  zi 

The  {q*}  are  due  to  impulsive  loads  wh^ch  are  known  in  the  blast  loading 

problem.  In  the  BR-1  code  {F},  [P},  [q*}  and  {q*}  are  known  at  time  t,. 

Hence,  {Aq*}.  and  fq*}.  can  be  determined  using  Eqs.  11-13.  For 
l+l  l+l 

our  situation,  {F}.  contains  {q*}  , which  is  unknown.  We  could 

1 i 

approximate  {q*}t  with  the  backward  finite  difference  form 

{4*}^  = (Aq*}t^  /At  (14) 

If  we  do,  then  {q*}  becomes  known  at  t , {F},  and  hence  {C},  can  be 

ti  i 

determined  at  t^,  and  the  procedure  used  in  BR-1  is  directly  applicable. 

On  the  other  hand,  if  we  express  fq*}*  in  the  central  finite  difference 

\ 

form 


. Consequently 


then  (F)  , and  hence  {C}.  , depends  upon  {Aq} 

*i  *i  *1+1 

{Aq*}  appears  on  both  the  left  and  right  hand  side  of  Eq.  11.  This 

Vl 

requires  a new  solution  procedure.  A detailed  study  of  the  accuracy 
and  numerical  stability  of  these  two  approacnes,  and  a third  approach, 
when  applied  tc  a single  degree  of  freedom,  damped  oscillator  is  pre- 
sented in  the  Appendix.  The  approach  where  {q*}  is  given  by  the  central 
difference  expression,  Eq.  15 » is  the  one  selected  based  upon  the 
accuracy  and  stability  properties  of  this  scheme.  Its  shown  in  the 
Appendix  that  the  maximum  value  of  At  for  a stable  solution  is  2/w,  where 
u,  is  the  highest  natural  undamped  frequency.  This  is  identical  to  the 
stability  limit  on  the  BR-1  procedure. 

Introducting  that  part  of  {F}  due  to  w given  by  Eq.  10a  into  Eq.  2 
results  in  the  modified  equations  of  motion 

CM]  {'4*}  + [D]  {4*}  = £C}  (16) 


Replacing  {q*}  with  the  conventional  central  difference  approximation 


-2  U*}t 

li 


(17a) 


is  equivalent  to  obtaining  {q*}t  from  Eq.  U with  {q*}  not  considered,  i.e. 


(17b) 


according  to  Eq.  12.  Substituting  Eq.  15  for  {q*}.  and  Eq.  17b  for 
{q*}.  into  Eq.  16  leads  to 


i+1 


(Aq*)t  = [M  + D (At/2)l-i  l [M  - D (At/2)]  {Aq*}t  + {C}t  I(l8) 


10 


* 


s 


which  is  equivalent  to  {Aq*}.  given  by  Eqs.  11  and  13  when  D is  a 

ti+l 

null  matrix  and  {q*}  is  not  considered. 

The  mass  matrix  [Ml  is  developed  in  BR-1  using  the  lumped  mass 
approach  and  is  given  by  (Eq.  A -97 » Ref.  1) 


[Mjl 


[Ml  = 


(19a) 


where  [Ml  is  a 6x6  matrix  given  by 


.,T 


Wr.  " £ W 1 


r nr  u n 


n a r 


(19b) 


and  [m  ] is  a diagonal  matrix  of  the  lumped  mass  at  the  a node  of  the 


nth  element.  Comparing  Eq.  10b  with  Eq.  19a  reveals  that  the  two 
matrices  [M  + D (At/2)]  and  [M  - D (At/2)]  occupy  the  same  nonzero  6x6 
locations  as  the  original  matrix  M.  Thus,  the  same  procedure  used  in 
BR-1  to  compute  [M]  1 can  be  used  to  compute  [M  + D (At/2)]  \ Its 
only  necessary  to  modify  the  elements  of  [M]y  by  the  addition  of  the 


damping  matrix  [D]  given  by  Eq.  10c.  The  other  necessary  change  is 
the  addition  of  the  Matrix  [M  - D (&t/2)l  as  a product  with  {Aq*}^  . 


Thus  since 


(20) 


Eq.  18  can  be  expressed  in  the  form 

(Aq*)  . = [M  + D (At/2)]-1  f[Mr  - Dr  (At/2)1  {Aq*}  (21) 

i+1  ' 1 

+ (^rtj  * r = 1»  2»*“R 
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C.  Program  Changes  and  Modification  Logic 


The  following  routines  of  the  IBM  version  of  BR-1  have  been  modified 
for  BR-1HR:  MAIN,  MEMBER,  MTERM,  OPIATE,  ST0RE,  DELTT , DEFDC,  REVIV1 
and  REVTV2.  Two  new  subroutines  were  created:  DFMASS  and  ADDAMP.  The 

core  size  was  increased  from  250k  bytes  to  290k  bytes. 

The  flow  of  the  logic  of  the  modifications  is  as  follows: 

1.  Compute  [M]  in  ST0RE  (no  change) 

2.  Compute  [D]y  in  ST0RE 

3.  Compute  [M]1  in  MTERM  (no  change) 

4.  Compute  maximum  time  interval  for  numerical  stability  in  DELTT 

based  on  [M]  (no  change) 

5.  Take  the  inverse  of  [M] to  gee  [M]y  in  DFMASS  using  INVS 

6.  Compute  [M  + Dr  (At/2)]  and  [M  - Dy  (At/2)]  in  DPMASS 

7.  Compute  [Mf  + Df  (At/2)]"1  in  FPMASS  using  INVS 

8.  Compute  [M  + Dr  (At/2)]"1  [Mf  - Dr  (At/2)l  in  DPMASS 

9.  Compute  [M  + D (At/A)]'1  r;j  - D (At/A)]  {Aq*}  in  ADDAMP 

r r r r 

10.  Compute  (Aq*}t  using  Eq.  21  in  DEFLX  (no  change) 

Lin 

The  phrase  "no  change"  means  that  the  original  procedure  was  used. 

When  no  damping  is  considered  the  modifications  and  additions  are 
bypassed. 
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III.  USER'S  INSTRUCTIONS  FCR  BR-1HR 


The  instructions  for  the  use  of  the  BR-1  code  are  given^  in  Ref.  2. 
All  of  the  Instructions  contained  in  that  volume  also  apply  to  the  modi- 
fied program  BR-1HR.  The  time  step  for  numerical  stability  of  BR-1HR 
is  identical  to  that  of  BR-1.  The  additional  instructions  required  to 
use  BR-1HR  are  as  follows: 


1. 


2. 


Problem  Control  Card  (page  4,  Ref.  2) 

IHR  (15,  Col.  66-70)  - IHR  = 0,  no  fluid  is  involved;  the 

original  BR-1  code  is  used.  IHR  » 1, 

(follows  IREV) 


fluid  is  involved,  the  modifications 
are  used. 

Rectangular  Panel  Card  (page  8,  Ref.  2) 

RH0CF  (E8.U,  Col.  55-62)  - RH0CF  is  the  product  of  Yf » the 

(between  BH0  and  Table  5^elflc  welSht’  and  c-  the 

sonic  velocity  of  the  fluid.  The 

units  of  VfCare  lbf  -in72  - sec?1 
If  the  panel  is  not  in  contact 

with  the  fluid,  RH0CF  = 0. 


1U 


/ 


7 


J 


s 


IV.  SAMPI£  PROBIEM 


A 8 imply  supported  square  plate  Is  subjected  to  a step  pressure  load 
of  the  form 


p = o 

p = P sin  — sin 

a a 

Due  to  symmetry,  only  one  quarter  of 
parameters  of  the  problem  are: 

E = 10.4  x 106  psi 
V = 0.0965  lb/ in? 
u * 1/3 
h * 0.1  in. 
a * 20  in. 

P = 0.01  lb/in? 


t * o 

t > o (22) 

the  plate  is  considered:  The 

- Young's  modulus 

- specific  weight  of  the  plate 

- Poisson's  ratio 

- thickness 

- length  and  width 


The  load  is  sufficiently  small  such  that  the  nonlinear  effects  are 
negligible.  The  plate  is  modeled  with  four  elements  as  shown  in  Fig.  1. 

The  equation  governing  the  damped  motion  of  the  plate  corresponding 
to  Eq.  16  is 

Dy\r  w + pew  = Psin  ^ sin  ^ (23) 

5 Aft 


where 


Eh' 


12(l-v2)  * 


+ 2 22 

ax  ax  ay 


JL 


.4 

+ 

ay 


and  g is  the  local  acceleration  due  to  gravity.  The  solution  to  Eq. 
23  is 
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E ■ 10.4  X I06p*i  »*l/3 

y«  0.0965  lb/in3  h * 0.1  In. 


FIGURE  I.  SAMPLE  PROBLEM 
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V - wit  [1  - e'Ctut  cos  (“\/l  - C2' (Ut  + J /cos  co}  (24) 

where 

tan  c = - C/V1  - C^1 

w = aS  sin  ex  ttjt 

st  7T4  a a 


when  the  plate  is  initially  at  rest.  The  displacement  at  the  center  of 
the  plate  given  by  Eq.  24  is  plotted  in  Figures  2,  3>  and  4 as  a function 
of  time  for  £ = 0,  0.666  and  270  corresponding  to  zero  damping,  less  than 
critical  damping  and  very  heavy  damping  respectively.  The  corresponding 


O 

values  of  gpc  for  the  fluid  are  0,  4,  and  1620  lbf  /(in  -sec).  Also 
plotted  in  Figs.  2-4  are  the  results  from  BR-1HR.  The  input  data  sheets 
and  the  print  of  the  input  data  are  given  in  Fig.  5*  The  execution  time 
on  the  IBM  360/67,  FORTRAN  IV  - Level  H,  was  8 min.  and  26  sec.  for  200 
time  steps  with  q = 270.  The  run  with  damping  not  considered  took 
essentially  the  same  length  of  time. 
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FIGURE  2.  TRANSVERSE  DISPLACEMENT  AT  NODE  9 VERSUS  TIME,  /■  0 (NO 


3. 


FIGURE  3. 


t (MSEC) 


TRANSVERSE  DISPLACEMENT  AT  NODE  9 VERSUS 
TIME,  ^ • 0 .666 


(s-0l 


1.4 


I.  2.  5 4.  5.  6 


* (MSEC) 


FIGURE  4.  TRANSVERSE  DISPLACEMENT  AT  NODE  9 VERSUS 
TIME  , J = 270 
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V.  SUMMARY  AND  CONCLUSIONS 


The  finite  element  digital  computer  code  BR-1  developed  by  the  Northrop 
Corporation  for  predicting  the  effects  of  internal  air  blast  on  typical 
combat  aircraft  skin- rib- stringer  structures  has  been  modified  to  include 
the  effect  of  compressible  fluid- structure  interaction.  The  fluid- structure 
Interaction  is  approximated  by  the  piston  theory  wherein  the  effect  of  the 
fluid  upon  the  structure  is  accounted  for  by  introducing  damping  to  the 
equations  of  motion  of  the  structure.  The  modified  code  is  called  BR-1HR. 
This  code,  in  conjunction  with  the  NWC  code  for  predicting  hydraulic  ram 
pressures,  can  be  used  to  predict  the  structural  response  of  aircraft  fuel 
tanks  subjected  to  penetrating  bullets  and  fragments. 

All  of  the  features  of  BR-1  exist  in  BR-1HR,  and  only  two  additional 
numbers  ere  required  for  the  input  data.  The  modified  code  is  operational 
on  the  IBM  360/67  in  FORTRAN  IV,  level  H,  and  requires  290K  bytes  of 
storage.  A sample  problem  was  executed  to  demonstrate  the  validity  of  the 
modified  code  for  zero  damping,  less  than  critical  damping,  and  very  heavy 
damping. 
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APPENDIX  - A STUDY  OF  THE  ACCURACY  AND  STABILITY  OF  SEVERAI  NUMERICAL 


INTEGRATION  SCHEMES  FOR  THE  TRANSIENT  RESPONSE  OF  HEAVILY  DAMPED 
STRUCTURES 

Many  studies  have  been  made  of  the  accuracy  and  stability  of 
numerical  integration  schemes  for  the  equations  of  motion  of  structural 
systems.  However,  most  of  these  studies  concentrate  on  the  response  of 
undamped,  or  lightly  damped,  systems.  Of  interest  here  is  the  response 
of  both  lightly  damped  and  heavily  damped  systems  since  both  kinds  of 
damping  can  occur  when  a structure  is  vibrating  in  contact  with  a 
compressible  fluid  (6). 

The  equations  of  motion  of  the  system  under  consideration  are  given 
in  matrix  form  by  Eq.  16.  Three  different  finite  difference  schemes  for 
the  numerical  solution  of  these  equations  are  considered  here.  Only 
explicit,  non-iterative  schemes  are  considered  due  to  the  fact  that  the 
BR-1  code  uses  an  explicit  solution  procedure.  Two  of  the  three  schemes 
are  based  upon  a two  variable  approach  UBing  {q*}  and  {v*} , where 

{q*}  = (v*)  (B-l) 

Thus,  Eq.  16  can  be  given  in  the  form 

{>}  = [M]-1  ({C}  - [D]  {v*})  (B-2) 

First  Scheme 

Substituting  the  first  order  approximations  for  (v*}  and  {q*} 

{v*}t  - (Wt  - {v*}t  ) /(dt)  (B-3a) 

zi  zi+l  zi 


(B-3b) 


{4*}ti+l  - ((<l*)ti+i  - {q*}ti)/(At)  - { AQ*}t ^ /fot) 
into  Eqs.  B-l  and  B-2  leads  to 

{v*}.  = [I  - M~^D  (At)]  (v*l.  + (At)  [M]"1  {Cl.  (B-Ua) 

*i+l  i 

(A4U  - At  {v»l.  (B-4b) 

ti+l  ti+l 

Eliminating  {v*l  from  Eqs.  B-4  results  in 

{A**}*  * [I  - M"1!)  (At)]  {Aq*}t  + (dt)2  [M]"1  (Cl  (B-5) 

i+1  i 

This  is  identical  to  the  scheme  used  in  BR-1  when  damping  is  not  considered. 
It  is  also  equivalent  to  the  scheme  where  the  acceleration  {q*]  is  approxi- 
mated by  the  conventional  central  difference  approximation,  Eq.  17a.  The 

two  variable  approach  given  by  Eqs.  B-4  may  be  more  desirable  than  Eq.  B-5 

2 

due  to  roundoff  error  considerations,  i.e.  (^t)  in  Eq.  B-5  is  a very 
small  number. 

Second  Scheme 

The  second  scheme  uses  Eq.  B-4a  and  the  simple  forward  Euler  approxi- 
mation for  {v*1  in  Eq.  B-l 

■ (At)  {v*}  (B-6) 

ti+l  *1 

in  place  of  the  backward  approximation  of  Eq.  B-4b.  The  solution  proce- 
dure is  to  compute  fv#l,  using  Eq.  B-4a  and  fAq*).  using  Eq.  B-6. 

Vl  1 ^i+l 

Third  Scheme 

The  third  scheme  uses  the  conventional  central  difference  approxi- 


mations  for  both  {q*1  and  {*4*1 , i.e.  Eqs  15  and  17a.  This  gives  Eq.  18, 
reproduced  here  for  convenience 

- [M  + D (At/2)]"1  ([M  - D (At/2)]  {Aq*^  + (At)2  (Cl^)  (B-7) 

is  also  identical  to  the  BR-1  scheme  when  damping  is  not 
A two  variable  version  of  this  scheme  is 

(AV*).  * [M  + D (At/2)]  ([M  - D (At/2)]  [av*]  + (a t)  {C}  ) (B-8a) 

li+l  ^i  li 

and 

{Aq*)t  = (At)  (Av*)+  (B-8b) 

li+l  *i+l 

2 

This  may  have  smaller  roundoff  error  than  Eq.  B-7  since  (^t)  has  been 
eliminated. 


This  scheme 
considered. 


The  Single  Degree  of  Freedom,  Damped  Oscillator 

The  equation  for  the  free  vibrations  of  the  single  degree  of  freedom, 
danced  oscillator  is 

mq  + dq  + hq  = 0 (B-9) 


Applying  the  three  schemes  described  above  to  Eq.  B-9  leads  to 


qt  - (2  - 2 « - fi  ) qt  + (1  - 2Cj)qt 


= 0 


i+1 


(B-lOa) 


1 0 

f ^ 

(At)v 

> - 

0 1 

q J 

► 

■ « 

L J 

*i+l 

-2 


1 - 2Ci 
1 1 


i-1 


r ^ 

(At)v 

► 

. q-t 

= 0 


(B-lOb) 


(1  + C(B)9+  “ ( 2 - (Jj  )q.  + (1  - C^q.  = 0 

i+1  i i-1 


(B-lOc) 
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where 


U)  * (a^)(d  tv  “ Vh/m  £ - d/  (2nvt0) 

according  to  Eqs.  B-5»  B-Ua  and  B-6,  and  B-7,  respectively. 
The  solution  to  Eq.  B-9  can  be  given  in  the  form 


(B-ll) 


where  A is  an  arbitrary  constant  and  the  superscript  i denotes  the  ith 
power.  The  solution  to  the  three  dif^-rence  schemes  for  an  arbitrary 
set  of  initial  conditions  can  be  obtained  by  assuming 


<L.  “ Ax1  (B-12a) 

^i 

(At)v.  -Bx1  (B-12b) 

*i 

where  x»  A and  B are  unknown  constants.  Substituting  Eqs.  B-12  into 
Eqs.  B-10  and  solving  for  \ for  each  scheme  lead  to 


1 - C5 


S2/2 


T <5 


\/(  C + <5- 


;/2  Y 


(B-13a) 


X2  - 1 - C5  = u\TF~-  1 (B»13b) 

X3  - U - £ /2  x £ i2/4  - l ) /(I  + c<5)  (B-13c) 

When  the  discriminant  in  Eqs.  B-ll  and  B-13  is  positive  the  solution 
consists  of  danped  motion  only.  When  it  is  negative  the  motion  is 
damped  and  oscillatory.  Thus,  a zero  discriminant  defines  the  limit 
of  the  oscillatory  behavior. 
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The  accuracy  of  the  three  numerical  schemes  can  be  evaluated  by 
comparing  \^,  Xg  and  Xj  with  \c  for  several  values  of  q and  £ . The 

values  of  the  ratios  of  the  numerical  solution  to  xc  are  given  in  Table 

B-l  for  £ « 0,  0.5»  5»  and  500  and  ft  = 0.1.  This  value  of  corresponds 
to  a time  step  of  1 ^ sec  when  ^ ■ 100,000  rad/sec  or  a time  step  of  1 
msec  when  ^ « 100  rad/sec,  etc.,  i.e.  the  solution  is  computed  ten  times 
over  the  time  interval  l/(0  or  20rr  times  over  the  undamped  natural  period 
2n/uc  The  closer  the  ratio  in  Table  B-l  is  to  one,  the  closer  the 
numerical  eigenvalue  is  to  the  correct  eigenvalue. 

The  numerical  stability  of  each  scheme  can  also  be  determined  from 
Eqs.  B-13.  When  |x|>l  the  numerical  solution  will  be  unstable.  The 
upper  limit  on  £ for  stability  can  be  determined  for  a given  value  of  Q 
by  equating  lx|  to  one  and  solving  for  ,5.  When  the  discriminant  is 
positive  \ * -1  is  the  limiting  value  and  the  negative  sign  in  the  * 

/ 2 5 

applies.  When  the  discriminant  is  negative  |x|  -/x  + y where  x and 
y are  the  real  and  imaginary  parts  respectfully.  The  results  for  the 
limiting  values  of  for  an  oscillatory  solution  and  for  numerical 
stability  are  given  in  Table  B-2. 


0.995004  T i O.947665  x i O.989949,  0.999495, 

0.099833  0.082275  0.371615  0.2062x10 


Scheme 

Oscillatory  limit 

Stability  limit 

X1 

£ " 2(1-C) 

i - 2(  \J<?  +l'-C) 

x2 

c - 1 

£ - 2c,  c 4 1 

£ - 2c-  2 v^T » 
C*l 

x3 

£ - 2 \/l  -C2 

£ - 2 

TAB  IE  B-2  Limits  on  £ for  an  Oscillatory  Solution  and  a Stable 


Solution 

Note  that  \2  is  unstable  for  any  non-zero  value  of  £ when  £*0  and  that 
2 is  the  maximum  limit  on  £ for  all  three  schemes.  Also  note  that  the 
limit  on  is  the  same  as  that  on  the  BR-1  routine,  even  with  damping, 
and  that  this  is  the  least  restrictive  scheme.  Consequently,  based 
upon  these  accuracy  and  stability  considerations  the  third  scheme  is 


selected. 


